import os, sys
import pickle
import Well
import math
from math import *
import WellMap
import random
from PIL import Image
import pylab
from pylab import *
from numpy import *
import KMeans
import pyXLWriter as xl
import copy

""" This module deals with the actual analysis of the experimental data.
It contains methods for well-map image analysis, clustering and data output. """

# Color tuple constant declartions for RGB schema
WHITE = (255, 255, 255)
BLACK = (0, 0, 0)
RED = (255, 0, 0)
GREEN = (0, 255, 0)

def unpickleFromFile(f):
    """ Loads and returns the pickled object found at f. """
    unPickleF = open(f, "r")
    return pickle.load(unPickleF)

def colourPixelsM(imgPath, avgMatrix, outColour, intensityThresholdScalar=2, out="coloredMPic.tif", bgColour=(0,0,0), ignoreSinglePixels=True):
    """ This method walks through an intensity image found at imgPath and
    colours all pixels that exceed the avgMatrix value times the intensityThresholdScalar.
    Pixels that do satisfy this criteria will be added as outColour pixels in new image. """
    print "Opening image: " + imgPath + " ..."
    try:
        img = Image.open(imgPath)
    except IOError:
        print "Problem opening img: " + imgPath + ". Skipping image."
        return
    print "Image opened."
    print "Analysing image..."  
    size = img.size;
    outImg = Image.new("RGB", size) 
  
    for x in range(size[0]):
        for y in range(size[1]):
            curPixel = img.getpixel((x,y));
            curPixelVal = int(curPixel)
            if curPixel >= intensityThresholdScalar*avgMatrix[x,y]:
                outImg.putpixel((x,y), outColour)
            else:
                outImg.putpixel((x,y), bgColour)
    
    # Walk through image dropping singly coloured pixels
    for x in range(size[0]):
        for y in range(size[1]):
            cP = outImg.getpixel((x,y))
            if not hasColouredNeighbor(outImg, (x,y), outColour):
                outImg.putpixel((x,y), bgColour)
            
    print "Finished. Writing image out to file..."
    outImg.save(out, "TIFF")
    print "Colored image written out to: " + out + "."
    return outImg

def buildBGMatrix(bgImages, writeOutMatrix=""):
    """ This method builds an average matrix for the images in the tuple
    bgImages. Note that we make the assumption that all images are of the
    same dimensionality. A[i,j] corresponds to the average pixel value
    at i,j over each of the images in the tuple; where A is the returned
    matrix."""
    # Open up each image file
    imgTuple = []
    for iStr in bgImages:
        try:
            cImg = Image.open(iStr)
            imgTuple.append(cImg)
        except IOError:
            print "error! couldn't open image: " + iStr
            
    img = imgTuple[0]
    avgMatrix = zeros(img.size, 'float')
    print "Building average background pixel intensity matrix..."
    for x in range(img.size[0]):
        for y in range(img.size[1]):
          pixelTotal = 0
          for curImage in imgTuple:
              pixelTotal += curImage.getpixel((x,y))
          avgMatrix[x,y] = float(pixelTotal)/float(len(bgImages))
    print "Matrix constructed. Returning matrix:"
    print avgMatrix
    print "\n"
    # if an output file was specified; pickle and write out the object
    if not writeOutMatrix=="":
        pickleFile = open(writeOutMatrix, "w")
        pickle.dump(avgMatrix, pickleFile)
    return avgMatrix

def pointFallsWithinBounds(p, bounds):
    """ Returns true iff p is within the square defined by bounds, which is 
    assumed to be two points. """
    p1 = bounds[0]
    p2 = bounds[1]
    #return p[0] >= p1[0] and p[1] >= p[1] and p[0] <= p2[0] and p[1] <= p2[1]
    return p[0] >= p1[0] and p[0] <= p2[0] and p[1] >= p1[1] and p[1] <= p2[1]

def imgToWellMap(img, bgMatrix, gThreshold, rThreshold, saveImgPath="wellTest.tif", picklePath="wellMap.pickle"):
    """ Converts an image to a well map. 
    
    Feed this an image where white pixels are assumed to be part of a 
    'well' and black is bg. This makes the assumption that wells are
    3x3 and drops anything <= 1pixel 
    
    Requisite parmaters:
    img --- this is the PIL image from which a well map will be constructed.
    
    Optional parameters:
    saveImgPath --- this is where an image corresponding to the constructed
                    well map will be saved.
    picklePath --- this is the path where the constructed well object will be 
                    written (pickled) to.
    """
    
    # Constants to be used in the matrix
    NOT_VISITED=0
    WELL=1
    NOT_WELL=2
    
    print "Analysing image for wells..."
    size = img.size
    wellImg = Image.new("RGB", size)
    
    # Create a fresh WellMap
    wMap = WellMap.wellMap(size, bgMatrix, gThreshold, rThreshold)
    
    alreadyVisited = zeros(img.size) 
    for x in range(size[0]):
        for y in range(size[1]):
            curPixel = img.getpixel((x,y))
            if alreadyVisited[x,y] == NOT_VISITED:    
                # Make sure this isn't a singleton pixel. If it
                # is, drop it.
                if curPixel == WHITE and hasColouredNeighbor(img, (x,y), WHITE):
                   upperLeftOfWell = getCoordsForUpperLeftOfWell(img, x,y)
                   xPrime = upperLeftOfWell[0]
                   yPrime = upperLeftOfWell[1]
                   # Note that xPrime and yPrime are unique. Let's 
                   # use them as indices. 
                   w = Well.well((xPrime, yPrime))
                   wMap.addWell(w)
                   for row in range(3):
                       for column in range(3):
                           try:
                               alreadyVisited[xPrime+row, yPrime+column] = WELL
                               wellImg.putpixel((xPrime+row, yPrime+column), WHITE)
                           except:
                               # We just ignore it if we go 'out of bounds'
                               continue
                else:
                   alreadyVisited[x,y] = NOT_WELL
                   wellImg.putpixel((x,y), (0,0,0))
    print "Well matrix constructed."
    wellImg.save(saveImgPath, "TIFF")
    print "Map written out as image to: " + saveImgPath
    pickleFile = open(picklePath, "w")
    pickle.dump(wMap, pickleFile)
    print "Map object pickled (saved) to " + picklePath
    print "Map built with " + str(wMap.numWells) + " wells."
    return wMap

def hasColouredNeighbor(img, pixelCoords, colour):
    """ Checks for neighboring coloured pixels.
    
    Returns true iff there is at least one neighbor of
    the pixel at pixelCoords that is has the colour specified
    by the colour parameter. 
    
    Requisite parameters:
    img --- PIL image object that will be looked at.
    pixelCoords --- coordinates for origin pixel (for which neighbors
                    will be considered). 
    colour --- will return true if one of the pixels adjacent to the
                pixel at pixelCoords is this colour.
    """
    for x in (-1, 0, 1):
        for y in (-1, 0, 1):
            if not x==y==0:
                try:
                    if img.getpixel((pixelCoords[0]+x, pixelCoords[1]+y))==colour:
                        return True
                except:
                    continue
    return False

def getCoordsForUpperLeftOfWell(img,x,y):
    """ Given a single coloured pixel coordinate specified by x,y, this
    finds the coordinates for the upper left hand of the well the pixel belongs to.
    
    The (idealized) assumption is that wells, at the pixel level, look like: 
    ppp
    ppp
    ppp
    And (x,y) is one of these pixels. This returns the upper left index
    tuple for the well.
    
    Requisite parameters:
    img --- the PIL image to be parsed.
    x --- x-coord of the pixel of interest
    y --- y-coord of the pixel of interest
    """
    # We consider the pixel (x,y) to be the upper left iff
    # the row above and column to the left are non white.
    # Checking column to the left.
    columnToTheLeftEmpty=False
    rowAboveEmpty=False
   
    columnToTheLeftEmpty = not isPixelWhite(img,(x-1, y-1)) and not isPixelWhite(img, (x-1, y)) and not isPixelWhite(img, (x-1, y+1))
    rowAboveEmpty = not isPixelWhite(img, (x-1, y-1)) and not isPixelWhite(img, (x, y-1)) and not isPixelWhite(img, (x+1, y-1))
    if columnToTheLeftEmpty and rowAboveEmpty:
        # This is the base case.
        return (x,y)
    else:
        recurseOn = None
        if not columnToTheLeftEmpty and not rowAboveEmpty:
            recurseOn=(x-1, y-1)
        elif not columnToTheLeftEmpty:
            recurseOn=(x-1, y)
        else:
            recurseOn=(x, y-1)
        return getCoordsForUpperLeftOfWell(img, recurseOn[0], recurseOn[1])
       
def getEuclidDistanceBetweenClusters(clustering): 
    """ Returns the Eluclidean distance between the centroids. This currently assumes only two clusters.
    @TODO extend to general case of n-clusters (perhaps summing the pairwise distances?) """
    return getEuclidDistance(clustering[0].centroid.coords, clustering[1].centroid.coords)
    
def getEuclidDistance(a, b):
    """ Euclidean distance between a and b is sqrt(sum((a[i]-b[i])^2) for all i) """
    ret = 0.0
    for i in range(len(a)):
        ret = ret+pow((a[i]-b[i]), 2)
    return math.sqrt(ret)

def isPixelWhite(img,c):
    """ Returns true iff pixel at coordinate c in image img is white. """
    try:
        p = img.getpixel(c)
        return p == WHITE
    except:
        return False

def genClusterDataTable(wells, wIntensitiesOverTime, wMap, excelOut, n=16, reclusterings=100):
    """ Returns a tuple; the clusters and the respective average red neighbor (in proximity n) count """
    print "Generating cluster data table.."
    print "Number of reclusterings (the best clustering will be returned): " + str(reclusterings)
    pList = []
    for w in range(len(wells)):
        # Map the points to the wells
        p = KMeans.Point(wIntensitiesOverTime[w], wells[w])
        pList.append(p)
    kClustering=None
    try:
        kClustering = KMeans.kmeans(pList, 2, .001)
        # 1/29/07 -- Recluster n times, minimizing the sum of squares error
        lowestError, tryCount = -1,0
        while tryCount < reclusterings:
            curClustering = KMeans.kmeans(pList, 2, .2)
            curSumOfSquares = KMeans.sumOfAllSquares(curClustering)
            #curDist = getEuclidDistanceBetweenClusters(curClustering)
            if curSumOfSquares < lowestError or lowestError == -1:
                kClustering = curClustering
                lowestError = curSumOfSquares
                print "Lowest Sum-of-Squares error thus far: " + str(lowestError)
            tryCount += 1
    except:
        print sys.exc_info()
        print "Ooops! Problem clustering.. Most likely there were no greens identified."
    #print len(kClustering)
    clAvgs = []
    if not kClustering == None:
        workbook  = xl.Writer(excelOut)
        worksheet = workbook.add_worksheet()
        worksheet.write([1,0], "Well Coordinate")
        worksheet.write([1,1], "Cluster Number")
        worksheet.write([1,2], "Red Neighbor Count")
        
        #print kClustering[0]
        print "Num Clusters: " + str(len(kClustering))
        print "In 0: " + str(len(kClustering[0].points))
        print "In 1: " + str(len(kClustering[1].points))
        #for well in wells:
        totesRCount, totesCount = 0,0
        clAvgsGreens = []
        #print range(len(kClustering))
        r = 2
        print "Neighborhood size is: " + str(n)
        for k in range(len(kClustering)):
            greenClustAvg = 0
            clustAvg = 0
            # For each cluster
            for wP in kClustering[k].points:
                # For each well in the cluster
                worksheet.write([r,0], "(" + str(wP.reference.x) + ", " + str(wP.reference.y) + ")")
                worksheet.write([r,1], str(k))
                rCount = wMap.countNumOfRedWells(wMap.getKNearestNeighbors(wP.reference, n))
                totesRCount+=rCount
                totesCount+=1
                clustAvg+=rCount
                gCount = wMap.countNumOfGreenWells(wMap.getKNearestNeighbors(wP.reference, n))
                greenClustAvg+=gCount
                worksheet.write([r,2], str(rCount))
                r+=1
            clAvgs.append(float(clustAvg)/float(len(kClustering[k].points)))
            clAvgsGreens.append(float(greenClustAvg)/float(len(kClustering[k].points)))
        avgOverClusters = float(totesRCount)/float(totesCount)
        worksheet.write([0,0], str(k+1) + " clusters generated. Average red neighbors over all clusters: " + str(avgOverClusters) + " Average number of red neighbors in " + str(n) + "-well neighborhood for clusters: " + str(clAvgs))
        workbook.close()
        print "Data file written out."
        for c in range(len(clAvgs)):
            print "Avg red neighbor count for cluster " + str(c) + ": " + str(clAvgs[c])
            print "Avg green neighbor count for cluster " + str(c) + ": " + str(clAvgsGreens[c])
        print "Done.\n"
    return [kClustering, clAvgs]

def writeAvgGraphOutForAllClusters(clustering, outDirPath, wMap, imgDir, avgRedNCount, iStart=None, iCutOff=None):    
    """ This method takes a clustering and dumpes the average graphs of intensities over time. """
    # @TODO This method needs to be cleaned up; there is some redundancy here and things ought to be split into methods
    plots=[]
    # 1/17/07
    normedPlots=[]
    maxI=-1
    for k in range(len(clustering)):
        # For-each cluster
        clusterWellLs = []
        print "Building well list for cluster " + str(k) + " ..."
        for wP in clustering[k].points:
            clusterWellLs.append(wP.reference)
        print "Done. Generating average graph for cluster " + str(k) + " ..."
        wellIntensitiesOverTime = getWellIntensitiesOverTime(imgDir, clusterWellLs, wMap, imgStart=iStart, imgCutOff=iCutOff)
        avgGraphPoints = genAvgGraph(wellIntensitiesOverTime)
        #avgGraphPoints = genAvgGraph(getWellIntensitiesOverTime(imgDir, clusterWellLs, wMap, imgStart=iStart, imgCutOff=iCutOff))
        # 1/17/07 -- plotting in normalized space
        normWellIntensitiesOverTime = map(lambda x:pylab.norm(x), wellIntensitiesOverTime)
        normedWellIntensities = copy.copy(wellIntensitiesOverTime)
        for k in range(len(normWellIntensitiesOverTime)):
            normedWellIntensities[k]=map(lambda x:x/normWellIntensitiesOverTime[k], wellIntensitiesOverTime[k])
        normedAvgGraphPoints = genAvgGraph(normedWellIntensities)
        normedPlots.append([range(len(normedAvgGraphPoints)), normedAvgGraphPoints])
        
        plots.append([range(len(avgGraphPoints)), avgGraphPoints])
        if max(avgGraphPoints)>maxI:
            maxI = max(avgGraphPoints)
    # 12/7/06 -- modified to give cluster graphs the same scale.
    k=0
    for p in plots:
        ind = arange(len(p[0])+2)
        # We want to insert a spacer at the start of the graph; hence we bump the range values.
        # e.g., [0,1,2,...] = [1,2,3,...]
        xCoords = [x+1 for x in p[0]]
        # Shift the xtick labels to reflect current time slice.
        xLabels = [""]
        xLabels.extend([str(x+iStart) for x in p[0]])
        xLabels.append("")
        plot(xCoords,p[1],"b-s")
        # Constrain axis coordinates
        axis([0,len(p[0])+1,0,int(maxI)+5])
        # 1/18/07 -- pretty plotting x axis
        xticks(ind, adjustAxisLabels(xLabels))
        xlabel("Time")
        ylabel("Intensity")
        savefig(outDirPath + "/cluster_" + str(k) + ".PNG")
        print "Done. Cluster " + str(k) + " graph saved to file."
        clf()
        k+=1
    
    # 1/17/07 -- plotting in normalized space. 
    # @TODO this should have been made a method (or at least refactored
    # to avoid pasting the whole thing here).
    k=0
    print "Plotting in normalized space.."
    for p in normedPlots:
        ind = arange(len(p[0])+2)
        # We want to insert a spacer at the start of the graph; hence we bump the range values.
        # e.g., [0,1,2,...] = [1,2,3,...]
        xCoords = [x+1 for x in p[0]]
        # Shift the xtick labels to reflect current time slice.
        xLabels = [""]
        xLabels.extend([str(x+iStart) for x in p[0]])
        xLabels.append("")
        plot(xCoords,p[1],"b-s")
        # Constrain axis coordinates
        axis([0,len(p[0])+1,0,1.1])
        # 1/18/07 -- pretty plotting x-axis
        xticks(ind, adjustAxisLabels(xLabels))
        xlabel("Time")
        ylabel("(Normalized) Intensity")
        savefig(outDirPath + "/cluster_" + str(k) + "_NORMALIZED_.PNG")
        print "Done. (Normalized) Cluster " + str(k) + " graph saved to file."
        clf()
        k+=1
    print "Done."
    
    # Now let's generate our scatter plot visualizing the clustering
    # We're going to assume only two clusters here.
    clusterSymbols=['s','o']
    colors = ['r','b']
    # For graph legend
    handles,legendLabels = [],[]
    for p in range(len(plots)):
        ind = arange(len(plots[p][0])+2)
        xCoords = [x+1 for x in plots[p][0]]
        xLabels = [""]
        xLabels.extend([str(x+iStart) for x in plots[p][0]])
        xLabels.append("")
        handles.append(scatter(xCoords,plots[p][1], marker=clusterSymbols[p], c=colors[p]))
        legendLabels.append("Sender Neighbor Count: " + str(avgRedNCount[p])[0:5])
        axis([0,len(plots[p][0])+1,0,int(maxI)+5])
    xlabel("Time")
    ylabel("Intensity")
    # 1/18/07 -- pretty-plotting x-axis
    xticks(ind, adjustAxisLabels(xLabels))
    legend(handles, legendLabels, 'upper left', shadow=True)
    savefig(outDirPath + "/clusters_scatter.PNG")
    clf()
    print "Cluster scatter plot saved to file."
    # 1/17/07 -- Adding scatter for normed vectors
    # @TODO this too should be split into it's own method; this and the above block
    # are redundant.
    handles,legendLabels = [],[]
    for p in range(len(normedPlots)):
        ind = arange(len(normedPlots[p][0])+2)
        xCoords = [x+1 for x in normedPlots[p][0]]
        xLabels = [""]
        xLabels.extend([str(x+iStart) for x in normedPlots[p][0]])
        xLabels.append("")
        handles.append(scatter(xCoords,normedPlots[p][1], marker=clusterSymbols[p], c=colors[p]))
        legendLabels.append("Sender Neighbor Count: " + str(avgRedNCount[p])[0:5])
        axis([0,len(normedPlots[p][0])+1,0,1.1])
    xlabel("Time")
    ylabel("(Normalized) Intensity")
    # 1/18/07 -- pretty-plotting x-axis labels
    xticks(ind, adjustAxisLabels(xLabels))
    legend(handles, legendLabels, 'upper left', shadow=True)
    savefig(outDirPath + "/clusters_scatter_NORMALIZED.PNG")
    clf()
    print "(Normalized) Cluster scatter plot saved to file."

def genAvgGraph(listOfPointSets):
    lsPrime = []
    # We assume all point sets of of equal length
    for p in range(len(listOfPointSets[0])):
        total = 0
        for ls in listOfPointSets:
            total+=ls[p]
        lsPrime.append(float(total)/float(len(listOfPointSets)))
    return lsPrime

def adjustAxisLabels(origLabelSet, maxLen=20, oneLabelPer=5):
    """ Helper method for splitting up long label sets (for pretty-plotting) """
    newSet = origLabelSet
    if len(origLabelSet)>maxLen:
        newSet = map(lambda x: "", copy.copy(origLabelSet))
        x=0
        while x*oneLabelPer < len(origLabelSet):
            newSet[x*oneLabelPer] = origLabelSet[x*oneLabelPer]
            x+=1
    return newSet

def threePlotTimeVsMaxIntensVsRedCountData(wMap, imgDir, lastImage, graphOut, n=32, out=None):
    """ This is for matlab style plotting (this file can be ported into matlab."""
    xPoints,yPoints,zPoints=[],[],[]
    imgList=os.listdir(imgDir)
    imgList.sort()
    imgSet = imgList[0:lastImage]
    yPoints = [x+1 for x in range(len(imgSet))]
    yPointsSet = []
    wSet = genKRandomWellCoords(wMap, 50, greenWellsOnly=True)
    t = 1
    #notRedWells = filter(lambda x: x not in wMap.redWells, gSet)
    for imagePath in imgSet:    
        print "Opening " + imgDir+"\\"+imagePath
        image = Image.open(imgDir+"\\"+imagePath) 
        #for w in notRedWells[0:100]:
        #for w in wSet:
        for w in wSet:
            xPoints.append(wMap.countNumOfRedWells(wMap.getKNearestNeighbors(w, n)))
            zPoints.append(wMap.getAvgIntensityForWell("green", w, image))
            yPointsSet.append(t)
        t+=1
    outF = open(out, 'w')
    outF.write("Red Neighbors\n")
    outF.write(str(xPoints))
    outF.write("Intensities\n")
    outF.write(str(zPoints))
    outF.write("Times\n")
    outF.write(str(yPointsSet))
    outF.close()
    
def plotMaxIntensVsRedCountData(wMap, imgDir, lastImage, graphOut, n=32):
    rCountLs=[]
    imgList=os.listdir(imgDir)
    imgList.sort()
    imgSet = imgList[0:lastImage]
    for imagePath in imgSet:
        print imgDir+"\\"+imagePath
        image = Image.open(imgDir+"\\"+imagePath)
        mostIntenseWells=findKMostIntenseWells(image, 25, wMap)
        #print mostIntenseWell
        rCount=0
        for k in mostIntenseWells:  
            rCount += wMap.countNumOfRedWells(wMap.getKNearestNeighbors(k[0], n))
        #rCount = wMap.countNumOfRedWells(wMap.getKNearestNeighbors(mostIntenseWell[0][0], n))
        print float(rCount)/float(25)
        rCountLs.append(float(rCount)/float(25))
    ind = arange(len(imgSet)+2)
    #xLabels = [str(x) for x in range(len(imgSet))]
    xCoords = [x+1 for x in range(len(imgSet))]
    axis([0,lastImage, 0, max(rCountLs)])
    ylabel("Sender Neighbor Count")
    xlabel("Time")
    title("Highest Responders vs. Sender Neighbor Count")
    plot(xCoords, rCountLs, 'r-s')
    #xticks(ind, xLabels)
    savefig(graphOut+"\\redcount_vs_mostintense_25.PNG")




"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
#
# 
#
# THESE METHODS NO LONGER USED
#
#
#
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""

def getAllColouredNeighbors(pixelLoc, img, bactPixelList, bactColour, alreadyVisited):
    x = pixelLoc[0]
    y = pixelLoc[1]
    neighboringPixels = ((x+1, y), (x-1, y), (x, y+1), (x, y-1), (x+1, y+1),
                            (x-1, y-1), (x-1, y+1), (x+1, y-1))
    for neighbor in neighboringPixels:
        # Bounds checking
        if neighbor in bactPixelList or neighbor[0]>=img.size[0] or neighbor[1]>=img.size[1] or neighbor[0]<0 or neighbor[1]<0:
            pass
        else:
            if img.getpixel(neighbor) == bactColour:
                # This neighbor is flouresced: e.g., the same bacteria
                # Append to list then recurse on neighbor
                bactPixelList.append(neighbor)
                alreadyVisited[neighbor[0], neighbor[1]]=1
                getAllColouredNeighbors(neighbor, img, bactPixelList, bactColour, alreadyVisited)
    return bactPixelList

def graphWellIntensitiesOverTime(imgDir, wellList, wMap, outDir="", writeGraphDataOut=True):
    """ For each well w in the provided wellList, this will look at each image in imgDir and 
    plot a graph showing the intensity of w over time. 
    
    Requisite parameters:
    imgDir -- This is a directory containing all the images of interest, which are assumed to be stored in sequential
                order over time. E.g., img1, img2 ... etc.
    wellList -- This is a list of well objects that give the coordinates for wells that will be plotted.
    wMap -- WellMap object corresponding to all the images in imgDir.
    
    Optional parameters:
    outDir -- The graphs will be written to this directory. If blank, this will try to write them
                to the imgDir + /Graphs directory, and will probably fail at doing so.
    writeGraphDataOut -- If true, as by default, this will write out the data points for each graph     
                (in tuple form: e.g., (time1, X)) for each well graph
    """
    # One graph for each well in wellList
    timeStamps = [] # This will be the x-axis of all the graphs
    wellIntensities = []
    for well in wellList:
        # Create an empty entry for each pixel
        wellIntensities.append([])
        
    if outDir == "":
        outDir = imgDir + "/Graphs"
    
    fList=os.listdir(imgDir)
    fList.sort()
    for image in fList:
        i = None
        print image
        try:
            i = Image.open(imgDir + "/" + image)
        except:
            print "Unable to open: " + imgDir + "/" + image
            
        if not i == None:
            for wIndex in range(len(wellList)):
                wellIntensities[wIndex].append(wMap.getAvgIntensityForWell("green", wellList[wIndex], i))
                
    # Now walk over each of the well lists
    w = 0 # Used for file naming only
    for well in wellIntensities:
        plot(range(len(well)), well)
        savefig(outDir + "/" + str(w) + ".PNG")
        w+=1
        clf() # Clear the image buffer
    if writeGraphDataOut:
        gFilePath = outDir + "/" + str(w) + ".graph.txt"
        gFile = open(gFilePath, "w")
        gFile.writelines(formatGraphData(wellIntensities, wellList))
        gFile.close()
    
    print "Finished writing out graph file!"

def getWellIntensitiesOverTime(imgDir, wellList, wMap, imgStart=None, imgCutOff=None):
    
    # One graph for each well in wellList
    timeStamps = [] # This will be the x-axis of all the graphs
    wellIntensities = []
    for well in wellList:
        # Create an empty entry for each pixel
        wellIntensities.append([])
    
    if imgCutOff == None and imgStart==None:
        fList=os.listdir(imgDir)
        fList.sort()
    elif imgStart==None:
        fList=os.listdir(imgDir)
        fList.sort()
        fList = fList[0:imgCutOff]
        print "Considering intensity from images 0 through " + str(imgCutOff)
    else:
        fList=os.listdir(imgDir)
        fList.sort()
        fList = fList[imgStart:imgCutOff]
        print "Considering intensity from images " + str(imgStart) + " through " + str(imgCutOff)
    for image in fList:
        i = None
        print image
        try:
            i = Image.open(imgDir + "/" + image)
        except:
            print "Unable to open: " + imgDir + "/" + image
            
        if not i == None:
            for wIndex in range(len(wellList)):
                wellIntensities[wIndex].append(wMap.getAvgIntensityForWell("green", wellList[wIndex], i))
    return wellIntensities

def colourPixels(imgPath, threshold, outColour, out="coloredPic.tif", bgColour=(0,0,0)):
    """ Colour pixels in the provided image and output the new image.
    
    This assumes the image at imgPath is a black and white (i.e., 
    not RGB image. The image will be written out to the file specified
    by the out parameter, with all pixels whose values exceed threshold
    being replaced by the outColor parameter; this should be a tuple (e.g.,
    (250, 0, 0) for all red. The bgColor by default is black, this can optionally
    be overwritten. 
    
    Requisite parameters:
    imgPath --- this is the path to the image that will be coloured.
    threshold --- pixels in the image at imgPath whose intensity is greater
                  than this will be coloured. 
    outColour --- pixels will be coloured this colour.
    
    Optional parameters:
    out --- path to output coloured image to.
    bgcolour --- non-coloured pixels will be coloured this.
    """    
                                                             
    print "Opening image " + imgPath  
    
    try:
        img = Image.open(imgPath)
    except IOError:
        print "problem opening img: " + imgPath + ". check path."
        return
    
    print "Scanning image..."    
    
    size = img.size;
    # We instantiate a new image because we need an RGB
    # instance.
    outImg = Image.new("RGB", size)
    for x in range(size[0]):
        for y in range(size[1]):
            curPixel = img.getpixel((x,y))
            curPixelVal = int(curPixel)
            if curPixel >= threshold:
                outImg.putpixel((x,y), outColour)
            else:
                outImg.putpixel((x,y), bgColour)
    outImg.save(out, "TIFF")
    print "Coloured image written out to: " + out + "."
    
def parseRGBImage(imgPath, out="pixelDataRGB.tif"):    
    # GFP is mapped to flourescent green pixels
    GFP = "02550"
    newGFP = (218,112,214)
    
    # Syto is mapped to red pixels
    Syto = "25500"
    newSyto = (255, 255, 0)
    
    GFPCount = SytoCount = 0
    GFPList = []
    SytoList = []
    
    print "Opening image " + imgPath  
    
    try:
        img = Image.open(imgPath)
    except IOError:
        print "problem opening img: " + imgPath + ". check path."
    outStr = ""
    print "Analysing image..."
    size = img.size;
    for x in range(size[0]):
        for y in range(size[1]):
            curPixel = img.getpixel((x,y));
            outStr+=" "
            curPixStr=""
            for k in curPixel:
                curPixStr+=str(k)
            if curPixStr == GFP:
                GFPList.append(("GFP"+str(GFPCount), x, y))
                GFPCount+=1
                img.putpixel((x,y), newGFP)
            elif curPixStr == Syto:
                img.putpixel((x,y), newSyto)
                SytoList.append(("Syto"+str(SytoCount), x, y))
                SytoCount+=1;
            else:
                outStr+=curPixStr
        outStr+="\n \n"
    print "Initial map analysis complete."
    print "Found: " + str(GFPCount) + " GFP Pixels and " + str(SytoCount) + " Syto pixels."
    print "Writing out coloured image to " + out + "..."
    img.save(out, "TIFF")
    print "Image has been written out."

def buildImgUsingWellMap(wellList, img, outPath="drawingWellsTest.tif"):
    """ Given a list of wells and an image, this creates a new image which
    colours *only* the well images. """
    newImg = Image.new("RGB", img.size)
    for w in wellList:
        try:
            # Let's just colour each well, ignoring other pixels in the
            # image. So for each well we colour the 3x3 neighborhood.
            for xPrime in (0,1,2):
                for yPrime in (0,1,2):
                    p = img.getpixel((w.x+xPrime, w.y+yPrime))
                    newImg.putpixel((w.x+xPrime, w.y+yPrime), p)
        except:
            # Some pixels might be out of bounds: we don't care.
            continue
    newImg.save(outPath, "TIFF")
    print "Coloured image written to " + outPath
    
def combineRedAndGreen(rImg, gImg, wells, outPath="", block=False):
    """ Interlaces red and green images. Will ignore any wells that
    contain both red and green pixels. "
    
    Walks through rImg and gImg examining all 3x3 wells specified
    in wells param. By default (blocks=false) it will decide if a given
    well is red or green; if a well is red it will copy over the pixels 
    from the input rImg to a new image - and if it is green it will
    do likewise wtih the green image. 
    
    Requisite parameters:
    rImg --- the image showing red bacteria.
    gImg --- the image showing green bacteria.
    wells --- list of well coordinates (each entry specifies an upper-left
                of a well.
    
    Optional parameters:
    outPath --- the interlaced image will be written out here if specified
    block --- all pixels in a well will be written out as either RED, BLACK or
                GREEN if this is true. Ignores intensities. 
    """
    
    # The images *should* be the same size...
    combImg = Image.new("RGB", max(rImg.size,gImg.size))
    for w in wells:
        if not block:
            putPixelsR = putPixelsG = zeros((3,3),object)
        wellIsRed = wellIsGreen = False
        rCount = gCount = 0
        try:
            # First we loop through to see what colour the well is
            for xPrime in (0,1,2):
                for yPrime in (0,1,2):
                    rPixel = rImg.getpixel((w.x+xPrime, w.y+yPrime))
                    gPixel = gImg.getpixel((w.x+xPrime, w.y+yPrime))
                    if rPixel != BLACK:
                        rCount+=1
                    #    if rCount > 1:
                        wellIsRed = True
                        if not block:
                            putPixelsR[xPrime][yPrime]=rImg.getpixel((w.x+xPrime, w.y+yPrime))
                    if gPixel != BLACK:
                        gCount+=1
                        #
                        # WE SHOULD DROP ALL SINGLE PIXELS BEFOREHAND
                        # Notice that we only count it as green if
                        # more than one pixel in the well is green.
                        #if gCount > 1:
                        wellIsGreen = True
                        if not block:
                            putPixelsG[xPrime][yPrime]=gImg.getpixel((w.x+xPrime, w.y+yPrime))
            # Now, if we've found a conflict - a well cannot be both
            # green and red - we just ignore it. 
            if (wellIsRed and wellIsGreen):
                # Ignore.
                colourWell = BLACK
            elif wellIsRed:
                colourWell = RED
            else:
                colourWell = GREEN
            # Now loop through again colouring the well.
            for xPrime in (0,1,2):
                for yPrime in (0,1,2):
                    # Even if we're not 'blocking' we want to 
                    # ignore 'wells' that are both red and green
                    if block:
                        combImg.putpixel((w.x+xPrime, w.y+yPrime), colourWell)
                    elif not colourWell==BLACK:
                        if wellIsRed:
                            combImg.putpixel((w.x+xPrime, w.y+yPrime), putPixelsR[xPrime][yPrime])
                        else:
                            combImg.putpixel((w.x+xPrime, w.y+yPrime), putPixelsG[xPrime][yPrime])
        except:
            continue
    if outPath != "":
        combImg.save(outPath, "TIFF")
        print "Combined image written to " + outPath

def colourWells(wells, imgSize, colour, outPath):
    print "Colouring a list of wells .."
    outImg = Image.new("RGB", imgSize)
    for w in wells:
        for xPrime in (0,1,2):
            for yPrime in (0,1,2):
                try:
                    outImg.putpixel((w.x+xPrime, w.y+yPrime), colour)
                except:
                    print "Pixel is out of bounds... ignoring"
                    continue
    outImg.save(outPath, "TIFF")
    print "Finished. Image output to " + outPath + "\n"
    
def getAllFeatureVectors(wellMap, greenDir, redImg, outDir=""):
    """ Builds all features vectors for all images in greenDir. This is meant to
    be used for a single experiment for which the WellMap object has
    been built (this is the wellMap param). """
    if outDir == "":
        outDir = greenDir + "/fVecOut"

    rIm = None
    try:
        rIm = Image.open(redImg)
    except:
        print "Not able to open red image at: " + redImg
        print "Returning."
        return None
    
    fList=os.listdir(greenDir)
    fList.sort()
    print greenDir
    print fList
    print len(fList)
    for gImage in fList:
        try:  
            print "Attempting to open green bacteria image at " + greenDir + "/" + gImage
            #Maybe I should colour them out for kicks..
            #colouredImg = colourPixelsM(dir + "/" + image, bgM, outColour, intensityThresholdScalar, outDir+"/"+image, bgColour)
            gIm = Image.open(greenDir + "/" + gImage)
            print len(wellMap.wellList)
            for well in wellMap.wellList:
                #print str(wellMap.numWells)
                fVec = wellMap.buildFeatureVector(rIm, gIm, wellMap.wellList[well])
                print fVec
                out = open(outDir + "/" + gImage + str(well) + ".pickle", "w")
                pickle.dump(fVec,out)
                print "Feature for " + str(w) + " built and pickled."
        except:
            continue
       
def pointFallsWithinBounds(p, bounds):
    """ Returns true iff p is within the square defined by bounds, which is 
    assumed to be two points. """
    p1 = bounds[0]
    p2 = bounds[1]
    #return p[0] >= p1[0] and p[1] >= p[1] and p[0] <= p2[0] and p[1] <= p2[1]
    return p[0] >= p1[0] and p[0] <= p2[0] and p[1] >= p1[1] and p[1] <= p2[1]
    
def graphRandomKWellIntensitiesOverTime(imgDir, k, wMap, outDir, boundBox=None):
    """ This will randomly choose k wells to plot over time. """
    randomWellList = genKRandomWellCoords(wMap, k, boundBox)
    graphWellIntensitiesOverTime(imgDir, randomWellList, wMap, outDir)
    
def graphRandomKWellIntensitiesOverTime(imgDir, k, wMap, outDir, boundBox=None):
    """ This will randomly choose k wells to plot over time. """
    randomWellList = genKRandomWellCoords(wMap, k, boundBox)
    graphWellIntensitiesOverTime(imgDir, randomWellList, wMap, outDir)

def genKRandomWellCoords(wMap, k, greenWellsOnly=False, redWellsOnly=False, boundBox=None):
    """ Generates k random well coordinates """
    randomWellList = []
    wellList = wMap.wellList
    while len(randomWellList) < k:
        well = None
        if greenWellsOnly:
            x = random.randint(0, len(wMap.greenWells))
            well = wMap.redWells[int(x)]
        elif redWellsOnly:
            x = random.randint(0, len(wMap.redWells))
            well = wMap.redWells[int(x)]
        else:
            x = random.randint(0, len(wellList))
            well = wellList[int(x)]
        
        if not boundBox == None: # If a boundbox is set
            if pointFallsWithinBounds((well.x, well.y), boundBox):
                randomWellList.append(well)
        else:
            randomWellList.append(well)
    return randomWellList


def formatGraphData(intensities, wellList):
    """ This will return a String pretty-printed of the graph points for the parametric well
        point list. Notice that each entry in intensities at index i is assumed to correspond to a list of
        length T containing time slice intensities for t in T for the well at index i."""
    s = " "
    # Write out horizontal 'header' giving well labels
    for w in wellList:
        s += "(" + str(w.x) + "," + str(w.y) + ") "
    s += "\n"
    print len(intensities[0])
    for i in range(len(intensities[0])): 
        s += str(i)
        for intensityList in intensities:
            try:
                s += " " + str(intensityList[i])
            except:
                print "ERROR!" 
                print str(i) + " out of range, here's the list: " 
                print intensityList
        s += "\n"
    return s
    
def findKMostIntenseWells(img, k, wMap):
    """ This will return a tuple of the k most intense wells in the parametric image. 

    Requisite parameters:
    img -- This is a (PIL) in which the k most intense wells will be found. Image must be gray scale.
    k -- Number of highest intensity wells in img to be returned.
    wMap -- The WellMap object corresponding to img.
    """
    # Record pixels as a tuple:
    #    ((x,y), intensity)
    # Where (x,y) is the pixel coordinate and intensity
    # gives the intensity level.
    kWells = []
    size = img.size
    for well in wMap.wellList:
        if len(kWells) < k:
            # If we haven't recorded any wells yet
            # then consider any arbitrary pixel
            # 'the brightest'. 
            # Notice that the average intensity function takes care of tNot regulation.
            kWells.append((well, wMap.getAvgIntensityForWell("green", well, img)))
        else:
            intensities = map(lambda x: x[1], kWells)
            minVal = min(intensities)
            wIntensity = wMap.getAvgIntensityForWell("green", well, img)
            if wIntensity > minVal:
                kWells[intensities.index(minVal)]=(well, wIntensity)
    return kWells



        
